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EXECUTIVE  SUMMARY 


This  report  is  a  summary  of  the  work  that  has  been  performed  on  the  implementation  of  global 
fracture  parameters  into  a  finite  element  (FE)  post-processing  computer  code.  The  work  is  a 
continuation  of  a  previous  study  involving  the  implementation  of  a  micromechanical  fracture 
model  into  the  LS-DYNA  user-defined  subroutines. 

Two  fracture  parameters  were  identified  and  implemented  into  the  post-processor;  the  strain 
energy  release  rate,  J ,  and  the  Weibull  stress,  cr^ .  The  calculation  of  these  parameters 

involved  parsing  the  output  data  of  the  selected  FE  code,  LS-DYNA,  including  element 
stresses,  strain  energies,  and  nodal  coordinates  and  kinematic  variables  for  an  arbitrary  number 
of  dump  times  in  a  specified  geometric  window.  The  calculation  of  J  involves  the  automatic 
determination  of  J-integral  contours,  which  was  coded  into  the  post-processor. 

The  calculation  routines  were  validated  using  a  linear  elastic  problem  for  the  J  -  integral 
calculator,  and  Beremin’s  [2]  data  for  the  Weibull  stresses.  The  calculated  results  were  in 
good  agreement  with  published  values.  Future  work  includes  the  continuing  testing  and 
improvement  of  the  post-processor. 
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1.0  INTRODUCTION 

This  report  is  a  continuation  of  previous  work  in  which  a  constitutive  model  for  steel  was 
developed  that  addressed  material  failure  in  both  the  brittle  and  ductile  regimes  [1].  For 
ductile  (high  temperature)  regions,  where  void  nucleation  and  coalescence  governs  failure,  the 
Johnson-Cook  plasticity  model  was  used  in  conjunction  with  Gurson  void  fraction  evolution. 
For  brittle  (low  temperature)  regions,  the  maximum  principal  stress  was  assumed  to  govern 
failure;  fracture  occurred  when  the  maximum  principal  stress  reached  a  critical  cleavage 
fracture  value. 

Although  a  critical  cleavage  fracture  stress  is  simple  in  concept,  it  is  not  widely  used  as  a 
practical  fracture  measure.  Investigators  have  tended  to  adopt  more  probabilistic  approaches 
such  as  the  Beremin  [2]  and  Master  Curve  [3]  methods.  In  these  methods,  a  scalar  fracture 
measure  for  cracked  regions  is  computed  or  experimentally  determined  or  both  for  a  number  of 
test  specimens.  This  fracture  measure  is  then  ranked  and  fit  to  a  cumulative  distribution  curve 
to  determine  the  probability  of  failure  at  a  specified  scalar  fracture  measure. 

In  this  project,  it  is  desired  to  adopt  the  probabilistic  methodology  to  compute  fracture 
initiation.  Two  measures  of  interest  include  the  strain  energy  release  rate,  J ,  and  Weibull 
stress,  .  The  original  tasks  were  as  follows: 

1)  Task  I:  Elastic-Plastic  Fracture  Mechanics  Parameter  and  Internal  Strain  Energy 
Density  -  a  methodology  to  compute  J  for  three-dimensional  finite  element  (FE) 
meshes  will  be  adopted  and  implemented  as  a  FE  post-processor.  This  will  also 
involve  calculation  of  the  internal  strain  energy  density,  W ,  in  the  EE  constitutive 
model  developed  in  the  previous  work.  The  J-integral  method  will  be  validated  using 
classical  test  problems  defined  by  the  Scientific  Authority  and  Martec. 

2)  Task  2:  Size  Dependent  Cleavage  Fracture  Modelling  -  the  calculation  of cr^  for 
specified  Weibull  parameters  will  be  added  to  the  FE  post-processor. 

3)  Task  3:  Application  of  the  Size  Dependent  DBT  Model  -  the  Beremin  model  will  be 
used  to  fit  Weibull  parameters  to  the  drop  tower  tests  analysed  in  [1].  The  J-integral 
method  may  be  used  to  rank  the  experiments  for  the  Beremin  model. 

After  the  completion  of  Tasks  1  and  2,  the  scope  of  Task  3  was  changed  slightly  because  of  the 
following  reasons: 

1)  It  was  unclear  how  to  apply  the  /-integral  and  Beremin  fracture  criteria  to  dynamic 
fracture.  Both  methods  are  more  suitable  for  determining  the  onset  of  brittle  fracture  in 
static  analysis. 

2)  Extensive  time  and  effort  had  been  spent  in  developing  a  crack  post-processor  for  ES- 
DYNA,  and  there  were  not  sufficient  funds  to  develop  a  methodology  for  dynamic 
fracture. 
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Task  3  was  rewritten  to  inelude  the  testing  of  the  Beremin  model  using  Beremin’s  original 
data,  using  speeimen  failure  loads  to  rank  the  experiments. 

In  order  to  eompute  these  global  fraeture  parameters,  a  erack  post-proeessor  for  LS-DYNA  has 
been  written  in  C++.  The  post-processor  reads  the  LS-DYNA  dSplot  file  for  the  desired  time 
window,  then  computes  either  Weibull  stresses  or  J-integrals  for  specific  times. 

Section  2  describes  the  methodology  used  for  the  3-D  J-integral  calculation  along  with  a 
validation  calculation.  Section  3  contains  the  Beremin  model  methodology  and  testing. 
Section  4  summarizes  the  work,  and  identifies  areas  of  future  study. 
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2.0  J  INTEGRAL  CALCULATION 

The  calculation  of  the  local  mechanical  energy  release  rate,  J ,  follows  the  procedure  outlined 
in  [4],  The  methodology  utilizes  a  domain  integral  approach,  where  an  average  value  of  J  is 
computed  for  the  3-D  domain.  The  calculation  is  performed  using  the  crack  post-processor 
described  in  Appendix  B. 


2,1  Preliminaries 

Figure  2.1  shows  the  crack  domain.  The  crack  coordinate  system  is  oriented  with  X2  normal 
to  the  crack  plane,  tangent  to  the  crack.  The  direction  X^  is  the  direction  of  crack 
extension.  The  calculation  volume  of  interest  is  denoted  as  V ,  and  is  enclosed  by  areas 
Aj  -  A5 .  The  crack  dimension  is  denoted  as  s  ;  V  extends  from  point  a  to  point  c  (s^  to  ). 
An  average  value  of  J  is  computed  at  point  b  . 


Figure  2,1:  J  Integral  Domain  Volume  and  Crack  Reference  Frame  (from  [4]) 

In  the  calculation,  the  following  is  assumed: 

1)  Body  forces  are  neglected  (other  than  inertial  forces). 

2)  Thermal  strains  are  neglected. 

3)  High  speed  crack  propagation  is  not  included. 

4)  Tensor  notation  is  used. 
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The  domain  integral  approach  utilizes  the  virtual  crack  extension  principle  [4],  In  this  method, 
the  contributing  terms  to  the  mechanical  release  rate  over  the  domain  from  a  to  c ,  J ,  are 
multiplied  by  an  arbitrary  weighting  function,  q  ,  as  shown  in  Figure  22. 


The  average  value  at  Z? ,  J ,  is  then  calculated  by  using  the  mean  value  theorem. 


J  = 


j  q{s)ds 


[2.1] 


Utilizing  assumptions  1-4,  ,  can  be  expressed  as 


^  I  1  I 

Ja-c  =  j  J{s)q{s)ds  =  j  qjPj,u^,-Wq,  --pu^  q,i+pUi  w,,i  q  dV ^  [2.2] 


where 

P  =  U'  Piola-Kirchoff  stress  tensor 

u,  u,u  =  displacement,  velocity,  and  acceleration,  respectively 
p  =  material  density 
W  =  strain  energy  density 
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and  the  commas  in  the  subscripts  imply  differentiation.  Since  the  1*'  Piola-Kirchoff  stresses 
are  used  here,  the  differentiation  is  with  respect  to  the  original  coordinates,  X  ,  rather  than  the 
current  coordinates,  x  .  The  strain  energy  density,  W  ,  is  computed  as 

W  =  j(T.jds,j  [2.3] 


where  cr  and  s  are  the  true  stress  and  true  strain,  respectively. 


2,2  Finite  Element  Formulation 

If  the  domain  is  discretized  into  a  number  of  finite  elements,  ,  Equation  2.2  becomes 


J 


a-c 


<=ir„V  ^  Ji 


[2.4] 


LS-DYNA  uses  one-point  integration  for  their  hexahedral  elements.  Using  this  integration 
rule  in  conjunction  with  the  standard  isoparametric  formulation.  Equation  2.4  is  now  written  as 


J„ 


=  Z  2 Mo  p{u){u]q,+  p{u)\u,}q 


[2.5] 


where  J  is  the  element  Jacobian  matrix.  The  element  coordinates,  kinematic  variables  and 
weighting  function  are  expressed  in  terms  of  trilinear  shape  functions,  N 


X,  =  (X(r,5,0){x,} 

[2.6] 

u  =  (^){u} 

[2.7] 

Ui  =  (X){u,} 

li,  =  (X){m,} 

[2.8] 

ii.  =  (x){ii',.} 

q  =  {N){q] 

[2.9] 

where  r,s,t  are  isoparametric  coordinates,  and  the  overbar  denotes  nodal  quantities. 
Appendix  A  contains  the  hexahedral  element  shape  functions. 

2,3  Element  Jacobian 

The  element  Jacobian  is  evaluated  using 
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[^] 


dX, 

dX^ 

0X3 

dr 

dr 

dr 

dX, 

dX^ 

0X3 

ds 

ds 

ds 

dX, 

dX, 

0X3 

dt 

dt 

dt 

Substituting  Equation  2.6  into  2.10  yields 


[^] 


S{M) ,  )  d{N) 

dr  ’  dr 

g('V)  (y  1  DjN) 

ds  '  ds 


d{N) 

dr 

d{N) 

ds 

d{N) 


dt 


{X,] 


2,4  Spatial  Derivatives 

Spatial  derivatives  are  computed  using  Equation  2.12 


'  0 

r  0 1 

0X1 

dr 

0 

0X2 

0 

■  =  [^r- 

0 

ds 

0 

[dX,\ 

Eor  the  weighting  function,  q  ,  Equation  2.12  yields 


dq 

0Xj 

dq 

0X, 

.  =  [j]-. 

dq 

t’l.i 

[0X3J 

Displacement  gradients  are  computed  similarly  to  Equation  2.13.  The 
derivatives  are  contained  in  Appendix  A. 


[2.10] 


[2.11] 


[2.12] 


[2.13] 


shape  function 
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2,5  1®^  Piola-Kirchoff  Stress  Tensor 

The  1*'  Piola-Kirchoff  stress  tensor  is  written  as 

[p]=|F|[alF]-"  [2.14] 

where 

[f]  =  deformation  gradient  tensor 
[cr]  =  Cauchy  stress  tensor 

Taking  the  transpose  of  both  sides  of  Equation  2.14,  and  noting  that  the  Cauchy  stress  tensor  is 
symmetric 

[Ff  =\F\[Fy[a]  [2.15] 

Equation  2.15  is  used  in  Equation  2.5. 


2,6  Deformation  Gradient  Tensor 
The  deformation  gradient  tensor  is  written  as 


5xj 

5xj 

axj 

aXj 

ax, 

ax. 

ax. 

ax. 

dX, 

ax. 

ax. 

dx^ 

ax. 

ax. 

aXi 

ax. 

ax. 

Substituting  the  transpose  of  Equation  2.13  into  Equation  2.16 


[2.16] 


[2.17] 
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2,7  Finite  Element  Domain  For  J-Integral  Calculation 
2. 7.1  Determination  of  Crack  Coordinate  System 

Figure  2.3  shows  the  crack  edge  of  interest,  defined  by  three  nodes.  The  normal  in  the  Xj 
direction,  n^,  is  previously  defined.  The  normal  parallel  to  the  crack,  is  taken  as  the 

direction  of  the  normalized  component  of  the  vector  between  node  1  and  node  2,  v  , 
perpendicular  to 


n,  = 


V  -  («  2  •  V  )v 

||v  -  («2  •  v)v 


[2.18] 


The  direction  of  crack  propagation,  n^,  is  then  taken  as  the  cross  product  between  and 


[2.19] 


1 


All  tensorial  variables  must  be  transformed  to  the  crack  reference  frame.  For  vectors  such  as 
coordinates,  displacements,  velocities  and  accelerations 

hi  =  kKjf}.  [2.20] 

and  for  stresses 

FI  =  [rlFlF]  [2.21] 

Where  the  subscripts  “c”  and  “g”  refer  to  crack  and  global  reference  frames,  respectively,  and 
[r]  is  a  transformation  matrix,  given  by 


[2.22] 
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2. 7.2  Determination  of  Weighting  Function,  q 

The  weighting  funetion,  q ,  ean  be  any  arbitrary  function  within  the  J-integral  domain,  and 
must  be  zero  on  the  domain  boundary.  An  easy  function  commonly  adopted  by  investigators 
is  to  set  g-  =  1.0  on  all  nodes  in  the  domain  interior,  and  q  =  0.0  on  the  domain  boundary. 


2,8  Determination  of  J-Integral  Contours 

Figure  2.4  shows  the  mesh  of  a  typical  crack  problem,  arranged  in  a  Cartesian  fashion.  The  J- 
integral  domain  is  always  two  elements  wide,  and  is  defined  by  a  number  of  contours  as  shown 
by  the  red  arrowed  lines  in  Figure  2.4.  Each  contour  consists  of  a  number  of  contiguous 
elements  that  encircle  the  previous  contours,  and  include  all  elements  contained  in  the  previous 
contours. 

Given  the  crack  tip  nodes,  the  J-integral  contours  are  determined  automatically  by  finding 
nodal  pairs  that  are  contained  in  the  adjacent  elements.  For  contour  1  in  Figure  2.4,  the  crack 
tip  nodes  are  used  as  nodal  pairs.  For  contour  2,  the  exterior  nodes  of  contour  1  are  used  as 
nodal  pairs  for  adjacency  checking,  and  so  on. 

For  this  implementation,  only  regular  Cartesian  arrangements  and  full  8  node  hexahedrons  are 
admitted.  Future  development  will  include  degenerated  elements  and  multiple  crack  tip  nodes. 


2,9  J-iNTEGRAL  CALCULATION 
The  calculation  method  proceeds  as  follows; 

1)  The  crack  tip  coordinate  system  is  defined  by  using  Section  2.7.1.  Coordinates, 
kinematic  variables,  and  stresses  are  transformed  to  the  crack  tip  coordinate  system. 

For  each  contour; 

2)  The  nodal  weighting  function,  q ,  is  applied  to  interior  and  boundary  nodes  using 
Section  2.7.2. 

3)  The  J-integral  is  evaluated  using  the  expressions  from  Sections  2.2-2. 6. 
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Figure  2,4:  Crack  Tip  Mesh  and  J-Integral  Contours 


2,10  J-iNTEGRAL  CALCULATION  EXAMPLE 

A  linear  elastic  example  problem  is  used  for  validation;  the  J-integral  is  computed  using  the 
crack  post-processor  described  in  Appendix  B.  In  this  example,  a  central  through-thickness 
crack  is  introduced  in  a  thin  steel  plate  of  finite  width  and  loaded  statically  perpendicular  to  the 
crack.  Figure  2.5  shows  the  crack  geometry.  The  strain  energy  release  rate,  J ,  is  computed  as 


J  = 


[2.23] 


where  K  is  the  crack  tip  stress  intensity  factor,  calculated  as  [7] 


K  =  PcrJKa 


na 

B  =  Jsec — 
\  W 


[2.24] 

[2.25] 


Figure  2.6  shows  the  finite  element  model  used.  An  elastic  modulus  and  Poisson’s  ratio  of 
200,000  MPa  and  0.3  were  used,  respectively.  The  plate  thickness  and  height  are  5  and  200 
mm,  respectively.  Ten  hexahedral  elements  were  used  through  the  thickness,  and  J  was 
calculated  at  the  plate  center.  An  end  load  of  1.0  MPa  was  applied. 
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G 


W  =  1 00  mm 
a  =  20  mm 


Figure  2,5:  Crack  Geometry 


Figure  2,6:  Finite  Element  Model 
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Table  2.1  shows  the  computed  J  for  each  of  the  contours.  The  average  J  for  contours  2-10 
was  3.9601x10''^  MPa-mm,  which  compares  favorably  (within  2%)  to  the  theoretical  value  of 
3.8833x10'^  MPa-mm. 


Table  2,1:  Calculated  J  Value  For  Contours 


Contour  Number 

J  (xlO"*  MPa-mm) 

%  Error 

1 

3.2622 

-15.99 

2 

3.8793 

-0.10 

3 

3.8693 

-0.36 

4 

4.0390 

+4.01 

5 

3.8997 

+0.42 

6 

4.0312 

+3.81 

7 

3.9761 

+2.39 

8 

3.9877 

+2.69 

9 

3.9679 

+2.18 

10 

3.9908 

+2.77 

Average  2-9 

3.9601 

+1.98 
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3.0  WEIBULL  STRESSES  AND  BEREMIN  MODEL  CALCULATION 

The  Beremin  model  [2]  is  a  methodology  that  eombines  finite  element  analysis  with 
experimental  results  in  order  to  estimate  a  measure  of  cleavage  fracture  stress  intensity  and  its 
uncertainty.  A  measure  of  cleavage  fracture  initiation  (such  as  displacement  at  fracture 
initiation  for  Charpy  V-notch  tests)  is  first  selected.  Normally,  tests  are  selected  in  the  brittle 
range  (low  temperature)  in  order  to  exclude  all  ductile  effects.  The  experimental  results  are 
listed  in  ascending  order  of  the  measure,  and  finite  element  analyses  are  performed  on  each 
test.  The  Weibull  stress  is  then  computed  and  fit  to  a  probability  distribution  curve. 

The  Weibull  stresses  are  computed  using  the  crack  post-processor  described  in  Appendix  B. 


3,1  Probability  Distribution  Curve 

The  failure  probability  of  any  given  test  can  be  obtained  using  a  two  parameter  probability 
distribution 


P.  =  \  -  exp 


f  \ 

G 

— 

V 

J 

y 

[3.1] 


where 


=  failure  probability 
=  Weibull  stress 

=  scaling  stress,  curve  fit  parameter  1 
m  =  Weibull  exponent,  curve  fit  parameter  2 


The  Weibull  stress  is  computed  using  a  summation  over  all  the  elements  in  the  crack  plastic 
zone.  Up 


cr,.,  = 


.Mm 


(  =  1 


V 

V^o  y 


[3.2] 


where 


Cj  =  maximum  principal  stress 
=  element  volume 

Fg  =  reference  volume,  can  be  taken  as  1 
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3,2  Probability  Function  Curve  Fit 


The  curve  fit  for  Equation  3.1  is  not  straightforward,  because  the  Weibull  stress  summation 
depends  on  the  m  parameter.  Usually,  m  is  assumed,  and  the  scaling  parameter,  cr„ ,  is 
calculated.  Equation  3.1  can  be  rearranged  to  form 


or 


where 


—  In 
m 


In 

V-^n) 


=  ln(CT„)- 111(0-.)  [3.3] 


y=x+h 


[3.4] 


y=— In 


m 


4^" 

V  V  s  J ) 


T=ln(cjJ 
Z)  =  ln(c7„) 


[3.5] 

[3.6] 

[3.7] 


3,3  Curve  Fitting  Procedure 

Initially,  test  data  is  ranked  according  to  a  monotonically  increasing  load  parameter  that  is  a 
direct  measure  of  cleavage  fracture  initiation.  Examples  of  load  parameters  that  are  commonly 
used  are:  load,  extension,  and  displacement  at  fracture  initiation.  Test  data  is  ranked  in 
ascending  order  of  the  load  parameter.  The  curve  fitting  method  then  proceeds  as  follows: 

1)  Erom  the  specimens,  a  failure  probability  is  assigned  to  each  specimen,  j ,  as 


^=Ur^,y  =  l,A'  [3.9] 

N 

2)  A  finite  element  analysis  is  performed  for  each  test,  and  the  Weibull  stress,  cr„, ,  is 
computed  using  Equation  3.2,  with  a  guessed  value  of  m  . 

3)  The  parameter  cr„  is  determined  from  the  x-intercept  of  curve  fitting  Equation  3.4. 
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3,4  Curve  Fitting  Example  -  Beremin  Data 


To  illustrate  the  procedure,  a  curve  fit  of  part  of  Beremin’s  [2]  data  is  shown  here.  A  multiple 
set  of  two  uniaxial  specimens  were  loaded  to  failure  at  low  temperature  (77  K).  To  check  the 
methodology  and  the  postprocessor,  a  set  of  finite  element  analyses  of  two  specimens, 
designated  AE0.2  Type  2  and  AE0.2  Type  3,  are  processed  to  obtain  the  Weibull  stresses  and 
curve  fits.  Eigure  3.1  shows  the  specimen  dimensions,  and  Eigure  3.2  the  EE  models  for  the 
Type  2  and  3  geometry.  One-quarter  symmetry  was  used  for  the  EE  models. 


Figure  3,1:  Specimen  Dimensions  for  AE0,2  Specimens  (obtained  from  [2]) 


Figure  3,2:  Specimen  FE  Models 
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Material  properties  were  obtained  from  Figure  3.3,  which  shows  the  temperature  dependence 
of  the  yield  and  flow  stress.  A  bilinear  stress-strain  law  was  assumed  for  the  FE  analyses.  The 
hardening  modulus  was  taken  as  the  difference  in  flow  stress  at  10%  strain  (e  in  Figure  3.3) 
and  the  yield  stress  divided  by  the  strain  difference.  The  data  shown  in  Figure  3.3  is  for  a 
different  heat  of  material  than  the  one  required;  Heat  A  is  required,  and  Heat  B  is  shown.  As  a 
result,  the  data  at  77°  K  is  scaled  by  the  yield  strength  difference  between  the  two  heats  (a 
reduction  factor  of  0.942).  The  ^MAT  PLASTIC  KINEMATIC  material  model  was  used; 
Table  3.1  shows  the  material  properties  used. 


TvmparMuraiKl 


Figure  3,3:  Yield  and  Flow  Stress  Dependence  on  Temperature  (from  [2]) 
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Table  3,1:  Material  Properties  For  Beremin  FE  Analysis 


Elastic  Modulus  (MPa) 

200,000 

Poisson’s  Ratio 

0.3 

Yield  Stress  (MPa) 

891.6 

Hardening  Modulus  (MPa) 

1449.0 

Each  FE  model  is  loaded  to  the  failure  loads  observed  in  the  experiments.  At  these  failure 
loads,  the  fraeture  post-processor  is  used  to  compute  the  Weibull  stress  assuming  an  m 
exponent  of  22.0.  Table  3.2  shows  the  sorted  Weibull  stresses  computed  from  the  FE 
analyses.  The  probability  of  failure  is  computed  using  Equation  3.9. 


Table  3,2:  Weibull  Stresses  -  Sorted  List 


Test  Number 

Type 

Weibull  Stress 

Pf 

1 

3 

1659.6 

0.02174 

2 

3 

1718.0 

0.06522 

3 

3 

1736.1 

0.10870 

4 

3 

1767.7 

0.15217 

5 

3 

1979.2 

0.19565 

6 

2 

2221.0 

0.23913 

7 

2 

2271.5 

0.28261 

8 

2 

2392.9 

0.32609 

9 

2 

2392.9 

0.36957 

10 

2 

2445.7 

0.41304 

11 

2 

2478.5 

0.45652 

12 

3 

2478.5 

0.50000 

13 

2 

2480.9 

0.54348 

14 

2 

2521.3 

0.58696 

15 

2 

2530.2 

0.63043 

16 

3 

2538.8 

0.67391 

17 

2 

2539.6 

0.71739 

18 

2 

2570.6 

0.76087 

19 

3 

2652.4 

0.80435 

20 

3 

2667.8 

0.84783 

21 

3 

2727.7 

0.89130 

22 

3 

2736.9 

0.93478 

23 

3 

2745.3 

0.97826 

Figure  3.4  shows  the  eumulative  plot,  along  with  the  curve  fit  from  Equation  3.3.  The  scaling 
parameter,  cr^,  was  determined  as  2517.7  MPa.  It  is  evident  that  the  curve  fit  is  less  than 
optimal.  There  could  be  several  reasons  for  this: 

1)  Selection  of  m  -  a  lower  value  would  rotate  the  central  portion  of  the  curve  counter¬ 
clockwise  and  fit  the  data  better  at  the  lower  tail.  This  can  be  seen  by  investigation  of 
Beremin’ s  original  data  in  Figure  3.5,  whieh  shows  essentially  the  same  behavior. 
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2)  Finite  element  model  -  in  the  crack  region,  it  was  found  that  relatively  few  elements 
resided  in  the  plastic  zone  at  the  lower  failure  loads.  A  more  refined  mesh  may  give 
more  consistent  results,  since  the  number  of  plastic  elements  would  increase.  It  is 
difficult  to  tell  what  refinement  Beremin  used,  since  the  mesh  data  is  not  available.  As 
well,  the  material  properties  used  in  the  Beremin  analyses  were  not  explicitly  specified. 
These  were  estimated  based  on  the  available  material  data. 


Figure  3,4:  Cumulative  Failure  Probability 


Nevertheless,  the  curve  fit  with  m  —  22  yielded  a  cr^  of  2518  MPa,  which  compared  favorably 
with  the  Beremin  computed  value  of  2570  MPa.  This  indicates  that  the  Weibull  stresses 
computed  using  the  post-processor  are  correct  values. 
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4.0  SUMMARY  AND  CONCLUSIONS 

A  post-processor  has  been  written  in  C++  code  to  compute  J-integrals  and  Weibull  stresses  for 
LS-DYNA  3-dimensional  hexahedral  elements.  The  post-processor  reads  LS-DYNA  binary 
fdes  (dSplot)  and  computes  J  and  cr^  for  an  arbitrary  number  of  contours  and  solution  dump 
times.  For  this  version,  the  crack  tip  is  assumed  to  consist  of  one  nodal  point  line,  and  the 
elements  are  8-noded  hexahedrons.  The  computational  modules  have  been  validated  on 
selected  test  problems. 

Although  great  progress  has  been  made  on  the  development  of  an  FE  global  fracture  post¬ 
processor,  it  remains  to  improve  and  test  the  software  for  larger  dynamic  problems.  The 
following  future  work  is  recommended. 

1)  Multiple  LS-DYNA  files  -  for  large  problems,  more  than  one  dSplot  file  is  generated. 
The  current  software  does  not  support  this. 

2)  Non-hexahedral  elements  -  incorporation  of  wedges,  tetrahedrons,  etc. 

3)  Tests  on  nonlinear  and  dynamic  problems  -  in  particular,  the  J-integral  calculation. 
This  can  be  performed  on  the  large  amount  of  data  collected  by  DRDC  Atlantic  for 
350WT  steel. 
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APPENDIX  A 

HEXAHEDRAL  SHAPE  FUNCTIONS  AND  DERIVATIVES 


APPENDIX  A:  ELEMENT  SHAPE  FUNCTIONS  AND  DERIVATIVES 


Hexahedral  element  shape  functions  can  be  expressed  as 

+  'I  ^)(l  +  +  h  «  =  1,8  [A.1] 

O 

Where  r,  s,  t  are  isoparametric  coordinates  that  extend  from  -1  to  +1,  and  i  is  the  node 
number.  Figure  A.l  shows  the  isoparametric  coordinate  and  node  numbering  convention 
utilized  by  LS-DYNA.  Table  A.l  shows  the  values  used  for  r.,  s,,  t.  according  to  node 
number. 


Figure  A.l:  Isoparametric  Coordinate  System 


i 

ri 
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-1 

-1 

-1 

2 
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-1 
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-1 

+1 

+1 

The  shape  function  derivatives  are  written  as 


dN, 


(l  +  .y.  s){l  +  1. 1) 


[A.2] 


dr  8 


[A.3] 


=  Y  (1  +  ^)(1  +  4  [A.4] 

or  8 


APPENDIX  B 

LS-DYNA  CRACK  POST-PROCESSOR 


APPENDIX  B:  LS-DYNA  CRACK  POST-PROCESSOR 


This  appendix  is  a  description  of  the  crack  post-processor  (hereafter  referred  to  as  CPP). 
CPP’s  function  is  to  parse  LS-DYNA  dSplot  files,  and  calculate  either  Weibull  stresses  or  J- 
integrals  for  a  selected  time  window.  The  program  requires  an  input  file  and  LS-DYNA 
dSplot  fide,  and  outputs  J-integral  values  and  Weibull  stresses.  The  input  file  is  in  keyword 
format  similar  to  that  used  by  LS-DYNA,  with  an  asterisk  (*)  preceding  any  keyword.  The 
keywords  can  be  in  any  order,  and  the  input  file  may  be  commented  with  the  first  character 
being  an  exclamation  point  (!).  The  input  is  free  format. 

B.l:  KEYWORD  DESCRIPTION 

*CALC 

Specify  the  solution  type 

Variable  Variable  Type  Comments 

calc  type  Integer  Calculation  type 


=  0  J-integral 
=  1  Weibull  stresses 


*CRACK 


Define  the  crack  orientation 


Variable 


Variable  Type 


Comments 


n  crack  nodes 


Integer 


Number  of  crack  nodes 
(3  for  now) 


crack_node[i] 
i=  1,3 


Integer 


Crack  node  numbers 


n2[i] 
i=  1,3 


Real 


Orientation  of  crack  n2  normal 
in  Figure  2.1 


*J  CONTOURS 


Specify  the  number  of  contours  for  J  calculations  and  properties 


Variable 


Variable  Type 


Comments 


contourtype 

ncontours 

partnum 


Integer 

Integer 

Integer 


Contour  type  -  not  presently  used 
Number  of  contours 
LS-DYNA  part  number  -  not  used 


*TIMES 


Specify  the  time  range  for  computation.  CPP  will  perform  calculations  for  dSplot  results 


within  this  time  window 

Variable 

Variable  Type 

Comments 

start  time 

Real 

Start  time 

end  time 

Real 

End  time 

*RHO 

Specify  the  material  density 

Variable 

Variable  Type 

Comments 

rho 

Real 

Material  density 

*MAT 

Specify  the  material  type.  For  a  material  type  of  zero  (elastic),  CPP  will  compute  the  elastic 
strain  energy  density  for  use  in  J-integral  calculations.  For  anything  else,  CPP  will  read  in  the 
strain  energy  density  from  the  dSplot  file 

Variable  Variable  Type  Comments 

mat  type  Integer  Material  type 

=  0  Linear  elastic 


^DYNAMICS 

Specify  the  dynamics  type.  If  zero  (static),  velocities  and  accelerations  will  not  be  included  in 
the  J-integral  calculations 

Variable  Variable  Type  Comments 

dyn  type  Integer  Dynamics  type 

=  0  Static 
=  1  Dynamic 


*WEIBULL  WINDOW 

Specify  the  geometric  window  that  Weibull  stresses  will  be  computed  over.  CPP  will  only 
include  elements  whose  centers  lie  within  the  window 


Variable 


Variable  Type 


Comments 


x_min[i]  Real 

i=  1,3 

x_max[i]  Real 

1=  1,3 

*WEIBULL  WINDOW 


Specify  the  parameters  for  Weibull  stress  calculation. 


Variable 

nWeibull 

m_exp[i] 

1=1,  n  Weibull 


Variable  Type 

Integer 

Real 


VO  Real 

B.2:  INPUT  FILE  COMPOSITION 


Window  minimum  coordinates 


Window  maximum  coordinates 


Comments 

Number  of  Weibull  exponents 
Weibull  exponents,  m 

Scaling  volume 


Table  B.l  shows  the  applicability  of  all  keywords  to  the  J-integral  and  Weibull  stress 
computations.  Any  keyword  marked  with  a  checkmark  must  be  included  in  the  selected 
computation  type. 


Table  B.l:  Keyword  Applicability  to  Selected  Analysis  Type 


Keyword 

J Integral 

Weibull  Stress 

*CALC 

V 

*CRACK 

V 

CONTOURS 

V 

*TIMES 

V 

*RHO 

*MAT 

^DYNAMICS 

*WEIBULL  WINDOW 

*WEIBUEL  PARAMS 

